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Abstract 

A  finite  difference  model  of  a  four-electrode  tissue  conductivity 
measurement  system  was  developed  and  shown  to  be  within  1 0% 
of  theory.  The  model  is  useful  for  explaining  the  behavior  of 
conductivity  measurement  electrodes  in  tissue. 

Keywords:  Conductivity  measurement,  current,  voltage,  electrodes 


Introduction 
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He  also  found  equations  for  an  infinite  half  space  of 
anisotropic  material,  longitudinal  (long)  and  transverse 
(trans)  to  the  fiber  direction: 
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Four  electrode  measurement  (4EM)  systems  have  been  used 
to  measure  conductivity  in  tissue  where  electrode 
polarization  would  cause  unacceptable  inaccuracy  in 
measurement  results  and  where  directional  measurement  of 
conductivity  is  required1 *  (1-3).  The  system  described  here  is 
intended  for  use  in  testing  conductivity  in  brain  tissues  in 
vitro. 

The  response  of  the  4EM  has  been  well  characterized 
in  homogeneous,  isotropic  media  as  (1,  3): 
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Steendijk’s  solutions  only  hold  for  specific  geometries.  A 
real  world  measurement  in  anisotropic  or  mixed  tissue  would 
necessarily  require  a  much  more  complex  solution,  if  such  a 
solution  is  even  possible. 
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Equation  (1)  assumes  the  four  electrodes  are  linearly 
arranged  with  two  high  impedance  potential  measurement 
electrodes  (V1,V2)  flanked  by  two  current  injection 
electrodes  (11,12).  It  also  assumes  the  electrodes  have 
uniform  spacing  (a)  and  are  surrounded  by  a  medium  with 
conductivity  a.  The  injected  current  is  from  a  point  source 
(II)  and  is  collected  by  a  point  sink  (12)  (Figure  1). 

In  addition  to  the  response  function  (Equation  1), 
Robillard  and  Poussart  (4)  also  characterized  the  field  of 
view  of  the  4EM  system  as  approximately  3  a  in 
homogeneous  isotropic  media.  In  other  words,  conductivity 
differences  in  media  outside  a  radius  of  3a  from  any  of  the 
electrodes  should  have  little  effect  on  the  response. 

Steendijk  (5)  found  solutions  for  4EM  systems  for 
several  conductivity  environments.  Starting  with  the  infinite 
space  of  isotropic  material  solution  (1)  (Equation  1),  he 
found  a  solution  for  an  infinite  half  space  with  electrodes 
placed  on  the  air/media  boundary: 


II  VI  V2  12 


Figure  1.  Electrode  configuration  of  a  4EM  system. 


1  Less  polarization  than  a  two-electrode  measurement  system  where  the  current  is  applied  by  the  same  electrodes  that  measures  the  potential 

and  directional  measurement  as  opposed  to  a  coaxial  measurement  system  that  measures  bulk  conductivity  but  gives  no  information  about 

direction. 
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This  paper  will  demonstrate  a  finite  difference  numerical 
solution  based  on  a  three  dimensional  matrix  of  conductivity 
tensors  to  support  any  combination  of  included  regions, 
limited  only  by  the  resolution  and  size  of  the  grid.  The 
tradeoff  between  either  a  finer  grid  size  or  increased  spatial 
resolution  is  processing  time. 


And  the  right  hand  side  of  Equation  (8)  is  simply  equivalent 
to  the  current  at  the  point  xyz,  Ixyz.  The  resulting  equation 
can  be  discretized: 


If: 


1 1 

1 1  ^  I  xyz 


(10) 


Materials  and  methods 

Dey  and  Morrison  (6)  developed  a  finite  difference  solution 
for  regions  of  differing  isotropic  conductivity  in  geological 
media.  This  paper  will  extend  their  solution  to  allow 
differing  anisotropic  regions.  The  solution  begins  with  a  few 
basic  equations: 


J  =  oE 


(4) 


E  =  — V4> 


(5) 


V-J  =  ^5  (xyz) 


(6) 


Equation  (4)  is  the  differential  form  of  Ohm’s  law  that  relates 
current  density  (J)  in  a  medium  of  conductivity  a  to  the 
electric  field  (E)  in  the  medium.  Equation  (5)  holds  because 
the  field  is  irrotational  so  E  is  the  gradient  of  the  potential 
(())),  and  equation  (6)  relates  the  time  derivative  of  a  point 
charge  (p)  at  location  xyz  to  the  divergence  of  J. 

Combining  these  equations  yields: 

r  n  dp  f 7 " 

—  V  '  \?'xyzV(Pxyz\  —  $ (xyz) 

where  the  quantities  oxyz  and  cj)xyz  are  conductivity  and 
potential  at  point  xyz. 

The  derivative  of  charge  with  respect  to  time  is 

equivalent  to  the  current  within  a  volume  element:  rr  =  —. 
H  at  av 

Integrating  both  sides  over  the  volume  element  (dV): 


To  simplify  the  equations,  we  utilize  a  uniform  spatial  grid 
throughout  the  matrix  where  the  voxels  have  the  same 
dimension  in  each  direction,  Ax  =  Ay  =  Az  =  h.  The  left- 
hand  side  of  Equation  (10)  is  integrated  by  breaking  into 
separate  equations  for  each  face  and  using  the  finite 
difference  approximation  for  the  normal  derivative. 

Dey  and  Morrison  estimated  the  conductivity  at  each 
voxel  face  by  averaging  the  conductivities  of  four  adjacent 
voxels.  In  order  to  account  for  anisotropic  conductivities,  we 
used  conductivity  of  the  adjacent  cell  in  a  direction  normal 
to  the  face  being  integrated  {onorm).  This  should  introduce 
minimal  error  as  long  as  the  change  in  conductivity  is  gradual 
local  to  the  voxel.  Care  should  be  taken  in  interpreting  the 
potentials  immediately  adjacent  to  a  boundary  between  two 
conductivities. 

The  discretized  solutions  for  each  face  are  summed  to 
comprise  the  surface  integral: 


(ii) 


where: 

//«dA=h2 

and  the  values  of  the  conductivity  (oadj)  and  potential  (4>adj) 
of  the  adjacent  voxel  are  defined  as  in  Figure  2.  Coupling 
constants  as  defined  by  Dey  and  Morrison  (6)  are  developed 
for  each  face  and  for  the  object  voxel,  refer  to  Figure  (2)  for 
definition  of  conductivity  for  each  face: 

C±  =  oi  *  h  (12) 


-  fff  V  ■  [crxyzV<pxyz\  dV  =  JJJ  ES(xyZ)dV  (8) 


The  left  hand  side  of  Equation  (8)  can  be  converted  to  a 
surface  integral  using  Green’s  theorem: 


-fff  V  •  [axyzV(J)xyz]  dV 

JJJ  AV 

I 
II: 


(9) 


=  —  IJ  aVcf)  ■  n  dA 
~ll‘.£dA 


where  dA  is  the  surface  area  element. 


C*2  —  @2  *  h 

C3  —  (J3  *  /l 

C4  =  04  *  h 

Q  —  &S  *  ^ 

C6  =  cr6*h 

Cself  =  —(fil  +  c2  +  C3  +  C4  +  C5  +  C6) 

All  the  analysis  to  this  point  applies  to  voxels  in  the  interior 
of  the  volume  of  a  conductive  region.  Because  the  size  of  the 
region  is  necessarily  limited,  some  attention  needs  to  be  paid 
to  the  behavior  at  the  limits  of  the  computational  volume. 
Voxels  at  the  boundary  of  the  volume  require  a  boundary 
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condition  to  allow  accurate  prediction  of  potentials  across 
the  boundary. 

The  simplest  boundary  condition  is  to  assume  the 
potential  at  the  boundaries  is  a  constant  (usually  zero),  as  if 
the  boundaries  are  at  infinite  distance  from  the  probe 
(Dirichlet  boundary  condition).  Implementation  of  the 
Dirichlet  boundary  condition  resulted  in  an  overprediction  of 
the  measured  potentials  by  25%  for  anisotropic  cases, 
presumably  due  to  interactions  of  fields  with  finite 
boundaries.  Although  Robillard  and  Poussart  (4)  found  a 
field  of  view  of  3 a  in  an  isotropic  medium,  an  anisotropic 
medium  may  alter  the  field  of  view  by  extending  the  range 
of  the  currents  in  the  low  conductivity  direction. 

A  mixed  boundary  condition: 


dc|) 

ccc|>  +  p  — 


f(x,y,z) 


(13) 


may  yield  a  more  accurate  prediction.  The  field  at  a  distance 
(r)  from  the  point  current  source  behaves  as  cf)  oc  i  so  the 
boundary  condition  becomes: 


d§  1 

—  oc  -  — er -n 
dn  r2 


(14) 


where  er  is  the  unit  vector  in  the  direction  from  the  point 
current  source  to  the  boundary  location  and  n  is  the  unit 
vector  normal  to  the  face  at  the  boundary,  e  ■  ?  =  rcos(0)  so 
the  boundary  condition  is: 


acf)  +  (3 


(J>cos(0) 

r 


f(x,y,z) 


(15) 


where  cos(0)  is  the  cosine  of  the  angle  between  the  vector 
from  the  point  current  source  to  the  boundary  voxel  and  the 
normal  direction  of  the  face  on  the  boundary.  And  r  is  the 
distance  from  the  source  to  the  voxel.  This  definition  greatly 
complicates  the  execution  of  the  model  because  creation  of 
a  new  finite  difference  matrix  is  necessary  for  each  source 
location.  To  simplify  model  creation,  we  can  assume  the 
source  is  always  at  the  middle  of  the  geometric  model 
volume.  This  causes  minimal  error  when  the  source  is  near 
the  center  of  the  volume  (6). 

As  a  result  of  the  mixed  boundary  condition,  the 
potential  for  faces  on  the  boundary  is  assumed  to  be  zero  (for 
instance,  Cx  =  0)  and  Cself  (Equation  1 6)  gets  an  extra  term. 
(Note:  Face  6  is  opposite  Face  1) . 

The  subscripts  ctr  and  bnd  refer  to  the  center  of  the 
geometric  model  volume  and  the  voxel  on  the  boundary, 
respectively.  Because  the  Cxbnd  terms  are  zero  except  for  the 
face  on  the  boundary,  the  inclusion  of  the  extra  term  only 
slightly  increases  the  complexity  of  the  model. 


CSeif  -  _(Ci  +  C2  +  C3  +  C4  +  C5  +  C6)  (16) 

“  (Cibnd  "f  C2bnd  4~  C3bnd 
"f  C4bnd  Csbnd  4“  Cebnd) 


Where 


Cibnd  —  C6  * 


_ lzctr  zbndl _ 

V(Xctr-Xbnd)2  +(yctr-ybnd)2  +  (zctr-zbr 


So  the  discretization  of  Equation  (10)  is: 


Ci^i  +  C2<t>2  +  C3c|)3  +  C4cf>4  +  C5c|)5  +  C6c|)6  (17) 

4"  Cself^fyiJX)  —  ^xyz 


This  equation  can  be  written  in  matrix  form  as: 

[CM]  =  [S]  (18) 

Where  [C]  is  the  matrix  of  coupling  coefficients,  [(()]  is  the 
vector  of  unknown  potentials  and  [S]  is  the  source  vector 
which  is  zero  everywhere  except  at  the  location  of  the  point 
current  sources. 


a  cube. 

In  order  to  construct  the  finite  difference  matrix,  the  3- 
dimensional  physical  matrix  must  be  mapped  to  a  2- 
dimensional  matrix  where  each  column  and  each  row 
represent  one  of  the  unknown  potentials.  For  instance,  a 
10x10x10  cell  region  (1000  cells  or  voxels)  would  map  to  a 
2-dimensional  matrix  with  1000x1000  elements  represent¬ 
ing  1000  equations  with  1000  unknowns  that  must  be  solved 
simultaneously  to  determine  the  potential  in  each  cell.  Each 
row  of  the  matrix  is  the  equation  for  one  unknown  potential 
and  each  column  element  represents  one  unknown  of  the 
equation. 
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Figure  3.  Top  right  had  comer  of  finite  difference 
matrix.  C  is  the  multiplier  of  the  unknown 
potential  at  each  adjacent  face. 

Figure  3  shows  the  top  right-hand  comer  of  the  finite 
difference  matrix  representing  the  system  of  n  equations  and 
n  unknowns  where  the  Cs  are  the  coupling  coefficients  with 
the  adjacent  cell.  Because  the  potential  in  a  cell  only  depends 
on  the  six  adjacent  potentials  and  its  own  potential,  the  finite 
difference  matrix  is  sparse  and  can  be  manipulated  using 
iterative  techniques  in  Matlab  (MathWorks  Inc.,  Natick, 
MA)  to  solve  for  the  potential  in  each  cell.  The  full  Matlab 
code  is  available  at  the  journal  website. 

Several  aspects  of  the  Matlab  code  deserve  mention. 
The  input  to  the  Matlab  code  is  a  geometric  model  consisting 
of  a  3  dimensional  array  of  diagonalized  conductivity 
tensors.  The  implementation  assumes  the  grid  spacing  to  be 
the  same  in  all  directions  (Ax  =  Ay  =  Az  =  h).  The  largest 
practical  size  of  the  grid  is  100^100x100  voxels  on  a  quad- 
core  MacBook  Pro.  With  a  grid  spacing  of  0.1  mm,  this 
provides  a  space  of  10  mm  x  10  mm  x  10  mm  that  is  large 
enough  to  contain  the  4EM  system  electrodes  with  a 
boundary  of  approximately  3  a.  The  grid  is  built  by  a  custom 
Matlab  function  (GM  builder.m)  using  a  mathematical 
description  of  the  medium.  The  medium  can  have  isotropic 
and  anisotropic  regions  and  could  be  imported  from  a 
diffusion  tensor  image  to  calculate  the  coupling  coefficients 
if  the  diffusion  tensor  is  assumed  to  be  proportional  to  the 
conductivity  tensor. 

The  most  computationally  expensive  part  of  the  code 
is  construction  of  the  finite  difference  matrix  (Matlab 
function  SetUpFDMatrix  impbound.m).  The  matrix  is  built 
one  line  at  a  time  using  logical  indexing  and  the  command 
sub2ind  that  converts  the  three  dimensional  coordinates 
(i,j,k)  to  a  serial  number  (1...N).  Matlab  references  matrix 
indices  in  row,  column,  slice  order  so  a  large  flat  matrix  can 
be  built  fairly  efficiently.  On  a  quad-core  MacBook  Pro  this 
step  can  take  from  20  minutes  to  one  hour  to  complete. 
Fortunately,  this  step  only  needs  to  be  accomplished  once  for 
each  geometry,  as  it  is  not  dependent  on  the  location  of  the 
sources  or  measurement  points. 

The  final  step  in  predicting  the  response  of  the  4EM 
system  is  to  specify  the  location  of  the  point  current  sources 
(one  positive  and  one  negative)  and  the  location  of  the 
measurement  points  (as  shown  in  Figure  1).  The  Matlab 
functions  SinglePoint.m  and  MultiPoints  angulareror.m 
compute  the  potential  for  a  single  probe  location  or  an  array 
of  probe  locations  within  the  medium.  The  functions  take  the 
center  of  the  probe  location  and  the  direction  of  the  probe 
axis  and  compute  the  location  of  each  current  point  source 
and  measurement  location  to  maintain  an  appropriate 
geometry.  In  a  physical  system  there  could  be  error  in  the 


positioning  of  the  probes,  so  the  function 
MultiPoints_angularerror.m  allows  specification  of  three 
angular  errors,  0X,  0y,  and0z.  The  function  can  show  the 
impact  of  having  the  probes  misaligned  with  the  interface 
directions  on  the  response  of  the  4EM. 

The  solution  of  the  finite  difference  equation  cannot  be 
accomplished  using  the  straightforward  Matlab  backslash 
operator  ([cf>]  =  [S]\[C])  due  to  memory  limitations  of  the 
computer.  Matlab  offers  several  methods  of  iterative 
solutions  for  sparse  matrices  that  make  the  solution  possible. 
The  General  Method  of  Residuals  (gmres)  (7)  is  a  robust 
method  that  is  successful  at  solving  the  system  regardless  of 
the  content  of  the  geometric  model.  However,  it  can  take  up 
to  an  hour  to  solve  the  equation.  The  Symmetric  LQ 
(symmLQ)  (8)  method  is  much  faster  for  some  geometries, 
but  fails  to  converge  if  the  geometries  are  asymmetric  such 
as  an  infinite  half-space.  The  method  we  selected  is  the 
Quasi-minimal  residual  (qmr)  (9)  method  because  it  is  very 
fast  (20-40  seconds)  and  converges  for  all  geometries 
attempted.  Note  that  all  the  methods  give  very  similar  results 
for  conditions  where  they  converge. 

Results 

The  finite  difference  model  can  be  compared  to  closed  form 
solutions  like  Steendijk  (5).  Equation  (1)  provides  a  solution 
for  an  infinite  space  with  isotropic  conductivity.  Equation  (4) 
provides  solutions  for  electrodes  on  the  surface  of  an 
anisotropic  halfspace.  Equation  (4)  can  be  modified  to 
predict  the  response  in  an  infinite  anisotropic  space  by 
multiplying  by  a  factor  of  1/2.  As  long  as  the  distance 
between  the  outer  electrodes  and  the  boundary  are  more  than 
approximately  three  times  the  electrode  spacing  (3xa)  the 
boundary  will  have  minimal  impact  on  the  measurement. 

Finite  difference  models  are  validated  based  on 
convergence  to  the  known  exact  value  for  the  solution  with 
decreasing  voxel  size  (h).  Figure  (4)  shows  a  graph  of 
predicted  potential  in  an  isotropic  infinite  space  with 
decreasing  6h’  as  compared  to  the  exact  solution  from 
Equation  (1).  Figures  (5)  and  (6)  show  graphs  of  predicted 
potential  versus  ‘h’  for  longitudinal  and  transverse  probe 
orientations  in  an  infinite  anisotropic  space  (o\  =  0.5,  at  = 
0.1).  These  plots  show  fairly  close  convergence  to  the 
analytical  solution  with  decreasing  voxel  size 
(Approximately  8%  difference  for  the  anisotropic  transverse 
cases,  approximately  2%  for  the  other  cases). 

Models  can  be  further  validated  by  comparing  them  to 
theory  over  a  range  of  values.  Figures  (7),  (8),  and  (9) 
compare  model  predictions  at  increasing  depth  with 
theoretical  predictions  (Equation  (4))  at  the  surface  of  an 
infinite  half  space  and  at  depth  (essentially  an  infinite  space). 
Note  that  the  finite  difference  model  slightly 
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Figure  4:  Convergence  plot  for  an  isotropic  infinite  space. 
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Figure  8:  Potential  on  the  surface  and  at  depth  in  an  anisotropic 
infinite  halfspace,  transverse  direction 
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Figure  5:  Convergence  plot  for  an  anisotropic  infinite  space, 
transverse  direction 


Figure  9:  Potential  on  the  surface  and  at  depth  in  an  anisotropic 
infinite  halfspace,  longitudinal  direction 
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Figure  6:  Convergence  diagram  for  an  anisotropic  infinite 
space,  longitudinal  direction 
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Figure  7:  Potential  on  the  surface  and  at  depth  in  an  isotropic 
infinite  halfspace 


underpredicts  (by  8%)  the  value  at  the  interface.  This  may  be 
inaccuracy  caused  by  the  simplistic  formulation  (kernel)  of 
the  finite  difference  model.  However,  at  depth,  the  model 
predicts  within  approximately  5%  of  the  theoretical  value. 

Accepting  that  the  finite  difference  model  is  a  fairly 
good  predictor  of  the  response  of  the  4EM,  a  number  of 
configurations  can  be  examined.  In  use,  the  4EM  is  advanced 
down  through  tissue  through  a  layer  that  is  thought  to  be 
anisotropic.  In  order  to  understand  the  response  as  the  probes 
traverse  the  anisotropic  layer,  a  geometry  model  was  built 
with  an  anisotropic  layer,  2mm  in  thickness,  flanked  by 
infinite  isotropic  layers.  Additionally,  an  isotropic  infinite 
space,  and  an  isotropic  infinite  halfspace  will  be  analyzed 
(Figure  (10)). 

Figure  (11)  shows  the  predicted  response  of  the  4EM 
system  as  the  probe  is  traversed  from  3  to  7mm  in  depth  (the 
z  direction)  in  the  center  (x  and  y  directions)  of  the  space. 
The  isotropic  media  in  these  analyses  have  an  electric 
conductivity  of  o  =  0.2  mS/cm  and  the  anisotropic  media 
have  conductivities  longitudinal  0\  =  0.5  mS/cm,  and 
transverse  at  =  0.1  mS/cm. 

The  response  of  the  4EM  system  probes  in  the 
anisotropic  layer  is  reduced  by  interactions  with  the 
surrounding  medium  as  compared  to  the  infinite  case.  In  the 
transverse  direction,  the  response  is  reduced  by  30.2%  and  in 
the  longitudinal  direction  the  response  is  reduced  by  20.2%. 
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Figure  1 1 :  Predicted  response  of  the  4EM  system  at  depths  from  3  to  7  mm 


Figure  12:  Probe  with  5  degree  positioning  error 


Figure  10:  Geometries  for  analysis  of  4EM  system  response.  Top 
figure:  infinite  space.  Central  figure:  infinite  halfspace  with 
probes  on  surface.  Bottom  figure:  layer  inclusion  in  infinite  space. 
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Figure  13:  Effect  of  5  degree  positioning  error  on  4EM  system 
response 
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Another  use  of  the  finite  difference  model  is  to  explore  the 
impact  of  probe  positioning  error  on  the  4EM  response. 
Figure  (12)  shows  the  probe  configuration  with  a  5 -degree 
positioning  error.  Note  that  the  probes  enter  the  layer  at 
different  probe  ‘z’  depths.  The  effect  of  this  can  clearly  be 
seen  in  the  graphs  in  Figure  (13).  The  error  causes  the 
response  of  the  4EM  system  to  be  reduced  by  4%  in  the 
longitudinal  direction  and  3.3%  in  the  transverse  direction. 

Discussion 

The  FD  model  converges  to  the  closed  form  result  and  agrees 
with  the  closed  form  solutions  across  a  range  of  values  within 
10%  in  all  cases.  Furthermore,  the  predicted  response  of  the 
4EM  system  is  reasonable  as  the  probes  traverse  down 
through  an  anisotropic  layer.  Future  work  on  this  effort  could 
include  measurement  of  anisotropic  media  of  known 
conductivity  to  further  validate  the  model. 

The  FD  model  will  be  extremely  useful  to  those  who 
need  to  directly  measure  conductivity  with  depth  in  tissue  if 
the  nature  of  the  tissue  structure  is  known  well  enough  to 
form  a  geometric  model.  It  is  also  useful  to  analyze  a  range 
of  tissue  and  geometric  properties  (e.g.  layer  thickness, 
anisotropy  ratio)  for  comparison  to  measured  data  and  to 
explain  4EM  system  response  in  anisotropic  heterogeneous 
media.  The  demonstrated  errors  caused  by  a  5 -degree  error 
in  probe  positioning  as  shown  in  Figure  13,  help  provide 
confidence  in  the  usefulness  of  the  4EM  system  in  a  practical 
measurement  situation. 
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